function C = distance_gpr(a, b)
[~, n] = size(a);

[~, m] = size(b);
 
mu = (m/(n+m))*mean(b,2) + (n/(n+m))*mean(a,2);
  
a = bsxfun(@minus,a,mu); b = bsxfun(@minus,b,mu);
  
C = bsxfun(@plus,sum(a.*a,1)',bsxfun(@minus,sum(b.*b,1),2*a'*b));

C = max(C,0);          
end
